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PACS 02 . 50 . Ey - Stochastic processes 
PACS 47 . 55 . Kf - Particle-laden flows 

PACS 77 . 84 . Nh - Liquids, emulsions, and suspensions; liquid crystals 

Abstract. - We investigate velocity probability distribution functions (PDF) of sheared hard- 
sphere suspensions. As observed in our Stokes flow simulations and explained by our single-particle 
theory, these PDFs can show pronounced deviations from a Maxwell-Boltzmann distribution. The 
PDFs are symmetric around zero velocity and show a Gaussian core and exponential tails over 
more than six orders of magnitude of probability. Following the excellent agreement of our theory 
and simulation data, we demonstrate that the distribution functions scale with the shear rate, the 
particle volume concentration, as well as the fluid viscosity. 
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Introduction. — To describe the statistics of complex 
systems, often probability distribution functions (PDF) 
arc utilized. These distributions have been found to be 
of non-Gaussian shape in numerous fields of physics, in- 
cluding astrophysics [1], flow in porous media [2], turbu- 
lence [3] , granular media [4-6] , or suspensions [7-9] . How- 
ever, the underlying processes are often not understood. 
In this letter, we focus on particularly important systems 
showing non-Gaussian velocity PDFs, namely sedimenting 
hard-sphere suspensions confined between sheared walls 
(see Fig. 1). They appear in river beds, blood examina- 
tions, industrial food production, the application of paint, 
and many more situations. Detailed experiments have 
been performed for more than a hundred years, but ques- 
tions about the microstructure or structural relaxations of 
the sediment are still not well understood. 



Numerous authors have found that the PDF of parti- 
cle velocities P(v) is not of similar shape as for an ideal 
gas, i.e., like a Maxwellian. Instead, P(v) can show a pro- 
nounced non-equilibrium shape, where the probability of 
high velocities is substantially larger [7,8]. In this letter 
we present a single-particle theory and simulations to show 
that such non-equilibrium distributions can be described 
as a consequence of an irreversible driving process, where 
particles on average gain energy by one mechanism, but 
loose energy by another one. Here, energy is gained from 
the shear or gravitational forces causing particles to col- 



lide. Contrarily, energy is dissipated due to viscous damp- 
ing. This causes P{v) to consist of a Gaussian core and 
exponential tails. Even though we focus on a well defined 
system here, the processes described are of general nature 
and can be applicable to ostensibly different setups. 

Experimentally, Rouyer et al. [10] studied quasi 2D 
hard-sphere suspensions and found a stretched exponential 
P(v) with concentration dependent exponents between 1 
and 2 corresponding to exponential distributions for high 
concentrations and Gaussians for small particle counts. 
These results contradict theoretical predictions of a tran- 
sition from exponential to Gaussian with increasing vol- 
ume concentration [7,8]. However, both experimental and 
theoretical studies do not present sufficient statistics over 
more than 2-4 decades. It is important to note that if 
one does not have enough data points for high quality 
PDFs, final answers on the nature of the function cannot 
be given. Indeed, in the process of analyzing our data we 
found that even stretched exponentials can fit PDFs with 
purely exponential tails and Gaussian centers if only two 
to four decades of probability are covered. But as soon 
as more data is added, the exponential nature of the tails 
becomes distinct and it is impossible to fit the whole PDF 
with a single function. 

Here, we overcome such limitations by presenting PDFs 
consisting of up to 10 10 particle displacements each - al- 
lowing a statistics superior to any previous work. Our 
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data does not show deviations from purely exponential 
tails over 6-8 decades of probability. We show that P(v z ) 
scales linearly with the shear rate, volume concentration 
and viscosity. 

The simulated system is a 3 dimensional setup as shown 
in Fig. 1. Top and bottom walls are at distance N z and 
sheared with shear rate 7 = 2v s \ lellI /N z . All other bound- 
aries are periodic. A body force / acting on the otherwise 
neutrally-buoyant particles can be added to mimic gravity. 
If turned on, this force causes strong density gradients in 
the system as can be observed in Fig. 1. We consider 384 
to 1728 initially randomly placed suspended particles of 
equal radius a corresponding to a particle volume concen- 
tration <f> between 6.8% and 30.7%. The particle Reynolds 
number Re = j2a/v is kept between 0.012 and 0.07. 
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Fig. 1: (Color online) Sketch of the simulation setup. 

Simulation method. — The simulation method is 
composed of a lattice Boltzmann solver (LB) for the fluid 
and a molecular dynamics (MD) algorithm for the motion 
of particles. This approach and recent improvements were 
originally introduced by Ladd and coworkers [11-14] and 
are well established in the literature [11-16]. Thus, we 
only shortly describe it here. 

The LB approach allows to calculate long-range hydro- 
dynamic interactions between particles, by utilizing a dis- 
cretized version of Boltzmann's equation [17]. Here, po- 
sitions x are discretized on a 3D lattice with 19 discrete 
velocities Cj pointing from a site to its neighbors. Ev- 
ery Ci is related to a single particle distribution function 
fi(x) which is streamed to neighboring sites at every time 
step. After streaming, a collision takes place where the 
individual fi(x) relax towards an equilibrium distribution 
/j 6q . Local mass and momentum density are given by mo- 
ments of the fi(x). The movement of suspended parti- 
cles is modeled by Newton's equation of motion and ap- 
propriate boundary conditions are imposed at solid/fluid 
interfaces to exchange momentum. We find a particle ra- 
dius of a=1.25 lattice sites sufficient since for larger radii 



P(v z ) does not change significantly anymore, while the 
computational effort increases substantially. Also, in low 
density simulations long-range hydrodynamic interactions 
dominate which are correctly reproduced even by small 
particles. For dense systems, exact lubrication forces be- 
tween particle pairs and between particles and walls are 
applied [11-13]. If many particles come close to each other 
(less than 0.1 lattice spacings), a cluster implicit method 
is used for updating forces in the MD algorithm [13]. This 
algorithm does improve the stability of the code and we 
carefully checked by comparing simulations with different 
volume concentrations and body forces that it has no in- 
fluence on the shape of the PDF. 

The simulation volume is 64a x 8a x 48a and 7 is varied 
between 2.3 • 10~ 4 and 1 ■ 10~ 3 (in lattice units). The 
fluid density is kept constant and the kinematic viscosity 
is set to v = 0.05 if not specified otherwise. A single 
simulation runs for 6.25 million LB steps, where during 
the last 5 million steps the z component of the velocity of 
every particle is gathered in a histogram to obtain P(v z ). 
All distributions are normalized such that J dv z P(v z ) = 1 
and Jdv z v^P(v z ) = 1 with the RMS velocity vf MS = 

Theoretical aspects. — Our theoretical model is 
based on the balance between viscous dissipation and 
shear forcing in steady state. The shear and the result- 
ing particle collisions are modeled by random, diffusive 
forcing. Due to this forcing, the velocity Vj of the jth 
particle changes as -jA- = £j where £j is a white noise, 
(£,•) = and (&(*)&(*')) = 2D6ij6(t - t'). In parallel, 
particles slow down because of the viscous fluid. In accor- 
dance with the traditional drag law, the velocity decreases 
according to dv/dt = —f3v, resulting in the exponential de- 
cay v(t) = v(0)e - ^' in the absence of forcing. We model 
this viscous damping by reducing each time unit At the 
velocity by a factor 77 = e~^At according to v — ► ryv. In 
a sheared fluid, there is a well defined time scale for re- 
encounters with the boundary, setting the time scale for 
the damping process. This damping process was used by 
van Zon et al. as a model for forced granular media [18,19]. 

The velocity distribution obeys the linear but non-local 
equation 



dP{v) _ D d 2 P(v) 



at 



d 2 v 



l - P ( V -)-P{v). 



The first term on the right hand side represents changes 
due to diffusive forcing and the next two terms represent 
changes due to viscous damping. In our theory, interac- 
tions between particles are represented through the ran- 
dom forcing process reflecting that a particle undergoes 
diffusion as influenced by all other particles. In steady 
state, the left hand side of Eq. 1 vanishes. We note that the 
shape of P(v) is independent of the diffusion constant D. 
Indeed, by making the scaling transformation v — > vj\J~D 
we can eliminate D and assume without loss of generality 
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that D = 1. The shape of the distribution depends on the 
dissipation parameter r\ alone. The moments 



M„ 



dv v n P(v) 



satisfy the recursion relation 

M„ = n(n-l)(l-77")- 1 M„_ 
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Since the distribution is symmetric. P(v) = P(—v), the 
odd moments vanish and starting with Mq = 1, the even 
moments are 



M 2 
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Of particular interest is the normalized 4th moment 

K = M4/M2 = 6/(1 + rj 2 ). (5) 

The distribution is close to exponential for strong drag, 
k> — ► 6 as Tj — ► and close to Gaussian for very weak 
drag, k — > 3 as rj — > 1. This is confirmed by the limiting 
behaviors of all moments 



2n 



(2n)l 



(2n- l)!!(l-r?)" n r)->l. 



(6) 



The leading large-velocity behavior can be derived by 
using a heuristic argument. For sufficiently large veloci- 
ties, the term i]~ 1 P (vij" 1 ) is negligible and hence, 



dv 



—P{v) = P(v) 



(7) 



Thus, P(v) has an exponential tail, P(v) ~ cxp ( — \v\). 
The prefactor can be obtained from the Fourier transform 
F(k) that equals an infinite series 

poo 00 

F{k) = / dve lkv P{v) = J[ [1 + fcVT 1 - (8) 

^° m=0 

This expression follows from the steady-state analog of 
(1), (1 + k 2 )F(k) = F(r/k). The simple poles at ±i closest 
to the origin imply an exponential decay, i.e., 



P(v)^A( v )eM-n), 



(9) 



when \v\ — * 00. Re-summation yields the residue to this 
pole, and in turn, the prefactor 



The exponential behavior is robust in the limit v — * 00. 
However, in the weak drag limit, 77 — > 1, the exponential 
behavior holds only for extremely large velocities. When 
r) — > 1, (77 = 1 — e, £ — > 0) we expand the denominator 



in A(rf). 
sum to 



Keeping only the dominant terms simplifies the 
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(11) 



and to leading order, the prefactor is A oc \ exp 
Therefore, under weak damping, there is a cross-over be- 
tween a Maxwellian behavior as follows from (6) and an 
exponential one, 




V<^£~ 



V~^>£~ 



(12) 



The two expressions match P(v) ~ exp ( — e _1 ) at the 
crossover velocity v ~ e" 1 . Interestingly, the crossover to 
a non-Maxwellian does not affect the leading behavior of 
the moments. 

In summary, the theory predicts the non-equilibrium 
shape of the PDF as an interplay between energy being 
injected by a diffusive thermostat and dissipation due to 
the fluid drag. In general, the high-velocity tail is expo- 
nential. The theoretical results shown later in this letter 
are given by a Monte Carlo solution of the steady state 
case of Eq. 1. In these simulations, N particles are char- 
acterized by a velocity Uj. The velocities change through 
two independent processes: damping and random forcing. 
In the damping process, the velocity is reduced by a fixed 
factor Vi — > rjVi. In the forcing process, the particle ve- 
locity changes by a random increment Vi — * vi + £ where 
£ has zero mean and a unit variance. The steady-state 
distributions were obtained using over 10 10 points from 
simulations with 10 8 particles. 




Fig. 2: (Color online) P(v z ) for / = 0.72 • 10 4 corresponding 



(10) to a Stokes velocity of v B — lA- 10 



= 13.6% and shear rates 
3.3-10 -4 , 6.7-10 -4 , and 1-10" 3 (iie=0.017, 0.034, 0.05). The 
solid line is the steady state solution of Eq. 1 (rj — 0.73). The 
lower inset shows the unsealed data, where higher 7 relate to 
wider P(v z ). The upper inset shows a linear fit with slope 0.21 
of vf MS (7) for the PDFs presented in the main figure as well 
as four additional datasets (symbols). 
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Results. — First, we consider suspensions with con- 
stant <p an d various shear rates under the influence of a 
body force /. The dependence of P(v z ) on 7 for three 
representative values is depicted in Fig. 2. P{v z ) is sym- 
metric and (v z ) = for all cases considered in this letter. 
As shown in the lower inset, the not normalized P(v z ) 
widen for higher 7. However, a very good scaling is ob- 
served: all normalized curves collapse onto a single one. 
In the upper inset we show the influence of 7 on w^ MS : as 
expected from the theory, 7 only sets a scale for the veloc- 
ity corresponding to a linear relation between 7 and ti^ MS . 
P{v x ) (after deducting the shear velocity) and P{v y ) also 
show Gaussian cores and exponential tails, although their 
width and height are different. Thus, the general shape of 
distributions in different directions is essentially identical 
and they are therefore not shown. To obtain an insight 
into the properties of P(v z ), we compute the cumulant 
k and find that for all simulation parameters studied in 
this letter it varies between 3.8 and 4.6. Knowing k, we 
can compute r] = y6/K — 1. Due to the large number of 
data points in our histograms, we calculate k for periods 
of 1 million time steps each and use the arithmetic average 
of the last 5 million time steps of a simulation run. We 
find that k varies by up to 10% within a single simulation 
which is of the same order as the difference of the individ- 
ual PDFs in Fig. 2. Thus, we average the different curves 
as well in order to obtain a value for the cumulant to be 
utilized for the Monte Carlo solution of the steady state 
case of Eq. 1. For the collapse in Fig. 2 we get 77 = 0.73. 
As depicted here, the solid line given by the theory and 
the simulation data excellently agree over the full range of 
six decades of probability. 

Next, we consider neutrally-buoyant suspended hard 
spheres (/ = 0)under shear. The shear rate is kept fixed 
and the particle concentration is varied between <f> = 6.8% 
and 30.7%. Due to hydrodynamic interactions, the par- 
ticles tend to move to the center of the system, i.e., to 
an area where the shear is low creating a depicted region 
close to the walls. However, this effect does not change 
the general shape of P(v). The corresponding normal- 
ized P{v z ) are presented in Fig. 3a. As depicted in the 
figure, all PDFs except for the lowest particle concentra- 
tion <j> = 6.8% (circles) collapse onto a single curve. At 
very low <j>, the tails of P{v z ) are still not fully converged 
due to the limited number of particle-particle interactions 
taking place within the simulation time frame. Again, the 
solid line in Fig. 3a is given by the steady state solution of 
Eq. 1 with 77 = 0.69 being obtained from the 4th moment 
of P(v z ). As before, simulation and theory agree very well. 
The full circles in Fig. 4 depict the dependence of vf MS 
on <j>. For concentrations of at least = 13.6%, vf- MS (cj)) 
can be fitted by a line with slope 1.8 • 10 -4 . The disagree- 
ment of the linear fit for low <j> is consistent with the not 
fully converged PDFs as shown in Fig. 3a. By keeping 
all parameters except the kinematic viscosity v constant, 
the dependence of v on P(v z ) can be studied. As demon- 
strated by the squares depicting the dependence of v^ MS 




Fig. 3: (Color online) P(v z ) for / = 0, 7 = 6.7 • 1(T 4 and 
<f) = 6.8%, 13.6%, 20.5%, 23.9%, and 27.3%, Re = 0.034 (a). In 
Fig. b), <j> is kept at 13.6% and the kinematic viscosity is set 
to v = 0.017, 0.05, and 0.1 (Re = 0.099, 0.034, 0.067). In both 
figures, all data sets collapse onto a single curve and the lines 
are given by the theory with r\ = 0.69. 

on v in Fig. 4, v^ MS and thus P(v z ) is independent of 
the viscosity. Thus, the steady state curve obtained for 
different volume concentrations is identical to the one for 
different v as shown in Fig. 3b. It would be interesting 

V 




Fig. 4: (Color online) vf MS in dependence of 4> (circles) and v 
(squares). Data corresponds to P(v z ) as in Fig. 3, but covers 
a wider range of <j> and v. Note the different a;-axes. 

to study the influence of the body force / on the shape of 
the PDF. However, / and the shear forces are in a sub- 
tle interplay since the height of the steady state sediment 
depends on both parameters and thus influences the local 
concentration. To investigate this behavior is beyond the 
scope of this letter. 

Conclusion. — To conclude, the non-equilibrium 
PDFs reported in this letter are a consequence of the irre- 
versible nature of the driving process. On average, parti- 
cles gain energy by external forces causing particle-particle 
collisions but lose energy by viscous damping. Unlike in an 
ideal gas, these two mechanisms are not interchangeable. 
In other words, one cannot reverse the arrow of time and 
observe the same behavior. Our theoretical model cap- 
tures this irreversibility through the competition between 
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two non-equivalent driving mechanisms: energy dissipa- 
tion through a multiplicative process and energy injection 
through an ordinary additive diffusive thermostat. 

The theory describes all aspects of the distribution as 
demonstrated by an excellent agreement with our coupled 
LB/MD simulations of sheared suspensions: the velocity 
distribution functions P(v z ) exhibit Gaussian cores and 
exponential tails over at least 6 orders of magnitude of 
probability. This finding is consistent with experimental 
results as given for example in [4]. We also note that the 
complete shape of the distribution can be characterized by 
a single parameter, the normalized 4th moment that has a 
one to one correspondence with the theoretical dissipation 
parameter. Further, we confirmed that P(v z ) scales lin- 
early with the particle volume concentration as well as the 
shear rate and is independent on the solvent's viscosity. 

While various authors report on transitions between 
Gaussian and (stretched) exponential tails [7,8,10], our 
results do not confirm such transitions. A likely reason 
for the previously published results might be the lack of 
sufficient statistics. If only two to four orders of magnitude 
of probability are covered, P{v z ) can be reasonably well 
fitted by pure Gaussians or distributions with (stretched) 
exponential tails as well. For a conclusive answer on the 
shape of the PDF, we found that at least five to six orders 
of magnitude of probability and a proper normalization are 
necessary. Experiments in strongly driven granular parti- 
cles suspended in a fluid provide evidence of exponential 
velocity distributions [4] . In order to close the question of 
the general validity of the findings presented in this letter, 
further high quality experimental or numerical data would 
be needed. 
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